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ABSTRACT 

We present a new version of the fully 3D photoionization and dust radiative transfer code, MOCASSIN, 
that uses a Monte Carlo approach for the transfer of radiation. The X-ray enabled MOCASSIN allows 
a fully geometry independent description of low-density gaseous environments strongly photoionized 
by a radiation field extending from radio to gamma rays. The code has been thoroughly benchmarked 
against other established codes routinely used in the literature, using simple plane parallel models 
designed to test performance under standard conditions. We show the results of our benchmarking 
exercise and discuss applicability and limitations of the new code, which should be of guidance for 
future astrophysical studies with MOCASSIN. 
Subject headings: radiative transfer — plasmas 



1. INTRODUCTION 

Photoionized environments characterize a wide range 
of astrophysical problems involving sources of X- 
radiation. With the advent of new technology used for 
instruments on board of (e.g.) XMM-Newton and Chan- 
dra, high resolution spectroscopy of such environments 
has become a reality. Paerels et al. (2000), for instance, 
observed the photoionized wind in Cygnus X-3 with 
Chandra High Energy Transmission Grating Spectrom- 
eter (HETGS) showing the discrete emission to be ex- 
cited by recombination in a tenuous X-ray-photoionized 
medium which is not symmetric with the source of the 
wind. Other examples include the detection of several 
recombination emission lines (from Fe xxvi at 1.78 A to 
N VI at 29.08 A), by Jimenez-Garate et al. (2005) in a 
50 ks observation of the bright X-ray binary Hercules X-1 
with Chandra HETGS. We also note the Chandra ACIS 
observations of the deeply eclipsing cataclysmic variable 
DQ Herculis by Mukai et al. (2003), who were able to 
pin down the origin of the soft X-rays from this system 
as being due to scattering of the unseen central X-ray 
source, probably in an accretion disk wind. 

A number of ID photoionization codes, including 
G. Ferland's cloudy (Ferland et al. 1998) and T. Kali- 
man's XSTAR (Kallman & McCray 1980, Kallman & 
Bautista 2001) continue to represent powerful analytical 
tools for the analysis of astrophysical spectra from the 
X-ray to the infrared regime. These codes are designed 
for diffuse, optically thin media, which may also be ir- 
radiated by a non-thermal X-ray continuum, and have 
been, for instance, applied to the modelling of Narrow 



Line Regions (NLRs) of Active Galactic Nuclei (AGNs). 
Several other codes have been developed which specialize 
in emission and reflection spectra from optically thick hot 
photoionized media, irradiated by a non-thermal contin- 
uum extending to the hard X-ray region, such as X-ray 
irradiated accretion disks (e.g. Ross & Fabian, 1993; 
Nayakshin et al. 2000; Dumont et al. 2000). 

To date, these and most other photoionization codes 
have numerically solved the equations of radiative trans- 
fer (RT) under the assumption of spherical symmetry or 
in plane parallel geometries, whereupon the problem is 
reduced to a ID calculation. While very few real X-ray 
sources are spherically symmetric, this approximation 
has been driven by the available computing power and 
the complexity of the multi-dimensional case. Mauche 
et al. (2004) developed a Monte Carlo code to investi- 
gate the radiation transfer of Lya, Hea, and recombina- 
tion continua photons of H- and He-like C, N, O and Ne 
produced in the atmosphere of a relativistic black hole 
accretion disk. This code, however, while accounting for 
Compton scattering and photoabsorption followed by re- 
combination, does not calculate the ionization state of 
the plasma. To our knowledge, there are currently no 
general, self-consistent and publicly available X-ray pho- 
toionization and dust RT codes capable of working in 
3D. Although the ID codes mentioned above are power- 
ful tools for the analysis of the pan-chromatic spectra of 
numerous astrophysical environments, their application 
is restricted to rather simplified geometries. 

The computational demand of realistic 3D simulations 
has recently come within the reach of low-cost clusters. 
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Taking advantage of this, the first self-consistent, 3D 
photoionization and dust RT code was developed for the 
IR-UV regime using Monte Carlo techniques (Ercolano 
et al. 2003a, 2005). The code, mocassin (MOnte CArlo 
Simulations of Ionized Nebulae) was designed to build 
realistic models of photoionized environments of arbi- 
trary geometry and density distributions, and can simul- 
taneously treat the dust RT. The code can also treat 
illumination from multiple point- or arbitrarily extended 
sources. The fully parallel version of MOCASSIN (docu- 
mented and publicly available) is well-tested for classical 
nebulae, according to standard photoionization bench- 
marks (Pequignot et al., 2001), and has been successfully 
applied to the modeling of H ii regions (e.g. Ercolano, 
Bastian & Stasihska, 2007) and planetary nebulae (e.g. 
Ercolano et al. 2003b,c, 2004; Gongalves et al., 2006; 
Schwarz & Monteiro, 2006; Wright, Ercolano & Barlow, 
2006). We now present a significantly improved version 
of the MOCASSIN code (version 3.00) which extends to 
the X-ray regime. In the tradition of previous releases of 
the MOCASSIN code, the X-ray version will also be made 
publicly available to the scientific community shortly af- 
ter publication of this article^ . 

In Section [2] we summarize the physical processes and 
atomic data added/modified in the new implementation 
and discuss the applicability and limitations of the code 
in its current form. In Section [3] we compare our code 
to established ID codes for benchmark tests, comprising 
emission line spectra from model NLR and from three 
thin low-density slab models illuminated by a hard con- 
tinuum. A brief summary is given in Section [J] 

2. THE X-RAY ENABLED MOCASSIN CODE 

2.1. Basic philosophy and underlying assumptions 

The original MOCASSIN code was developed in or- 
der to provide a 3D modelling tool capable of dealing 
with asymmetric and density and/or chemically inho- 
mogeneous media, as well as, if required, multiple, non- 
centrally located, point and extended sources of excit- 
ing radiation. The code can self-consistently treat the 
transfer of ionizing and non-ionizing radiation through 
a medium composed of gas and/or dust. The numeri- 
cal techniques employed and the physical processes con- 
sidered by the code are described in detail by Ercolano 
et al. (2003a, 2005). In brief, MOCASSIN locally simu- 
lates the processes of absorption, re-emission and scat- 
tering of photons as they diffuse through a gaseous/dusty 
medium. The radiation field is expressed in terms of en- 
ergy packets which are the calculation 'quanta'. The 
energy packets are created at the illuminating source(s), 
which may be placed anywhere inside the grid and be of 
point-like or diffuse nature. Their trajectories through 
the nebula are computed, as the packets suffer scatter- 
ings, absorptions and re-emissions, according to the lo- 
cal gas and dust opacities and emissivities. The packets' 
trajectories yield a measure of the local radiation field, 
from which the local photoionization and recombination 
rates, as well as the heating and cooling integrals, are 
determined. 

It is worth reminding the reader of the physical as- 
sumptions that both the previous and current versions 

^ Previous versions are available upon request from BE. 



of the MOCASSIN code rely upon, which also define its 
applicability and limitations. 

• All physical processes affecting the gas are in 
steady state. This implies that the atomic physics 
and heating and cooling timescales are short com- 
pared to those of gas-motion or to the rate of 
change of the ionizing field. 

• In some cases large optical depths may occur in 
the core of emission lines of astrophysically abun- 
dant ions. As described by Ercolano et al. (2005), 
resonant scattering is accounted for via an escape 
probability method. The underlying assumption is 
that photons in the line centre will be scattered 
many times in a region close to where they were 
originally emitted, until they finally escape through 
the line wings or are destroyed through continuum 
photo-absorption. While escape probability meth- 
ods are often used in photoionization models, it 
is well known that these may lower the accuracy 
of calculations in high density environments (e.g. 
Avrett & Loeser, 1988; Dumont et al. 2003). Fur- 
thermore, our current scheme assumes a static gas 
distribution, which would need to be modified to 
treat winds or accretion disks with very large shear. 

• The contribution from recombination is included 
for H-like, He-like and the lower Fe xxiv lines. 
Recombination, however, can dominate for several 
other lines in some extreme cases; in general, this 
process strengthens the 3s Fe L-lines compared to 
the 3d (Liedahl et al., 1990). Nevertheless, in 
the most commonly encountered cases, the contri- 
butions from recombination of He-like and H-like 
species are strongest. 

• Very high density environments are problematic 
and solutions may carry larger uncertainties. Apart 
from the line transfer issue mentioned above, a 
higher density limit of ^ 10^^ cm~^, at temper- 
atures of ^ 10** K, is imposed by our approxi- 
mate treatment of three-body recombination and 
coUisional ionization, which become important at 
higher densities (see also Ferland 2006). However, 
for highly ionized, hotter gas, e.g. ^ 10^ K gas 
in the corona of an X-ray binary, the high density 
limit is increased to ^ 10^^ cm~^. 

• A treatment for unresolved transition arrays 
(UTAs), observed in several spectra of AGNs (e.g. 
Sako at al. 2001), is not currently included in 
MOCASSIN. Calculations were presented by Be- 
har, Sako and Kahn (2001) for inner shell n=2-3 
photoexcitation of the 16 iron charge states Fe I 
through Fe xvi and their data will be included in 
a future version. 

2.2. X-ray extension 
2.2.1. Atomic Data 

mocassin's atomic database was updated in order to 
include the latest data releases. Details of mocassin's 
complete database for the original version are given in 
Appendix A of Ercolano et al. (2003a). Data updates 
(additions and replacements) implemented for version 
3.00 are as follows: 
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Free-bound emission for H i, He ii and He i uses 
the data recently calculated by Ercolano & Storey 
(2006). 

The radiative and dielectronic rates of Badnell 
et al. (2003), Badnell (2G06a,b), Zatsarinny et 
al. (2003, 2004a,b, 2006), Colgan et al. (2003, 
2004), Altun et al. (2004, 2006), Mitnik & Bad- 
nell (2004) have been included and are used as de- 
fault. The data cover all elements up to Zn in- 
cluding sequences up to Na-like electron target, all 
other species are treated with the previously avail- 
able data, except for species included in the Nahar 
(1997,1999) data set, if chosen by the user (see fol- 
lowing item). 

The total (radiative and dielectronic) rates of Na- 
har (1997, 1999) for recombination to all ions of 
carbon, nitrogen and oxygen and of Nahar (2000) 
for recombination to Si i, S ill, Ar V, Ca Vii and 
Fe XII can now be included, if explicitly chosen. 
As well as the recombination data set mentioned 
as the default in the previous item, all data used 
in previous versions of MOCASSIN are still available 
for use and compatible with the X-ray version. 

The fits of Verner (1996) to Opacity Project 
(Seaton et al. 1993) data for the photoionization 
cross-sections from all shells are used (the previ- 
ous versions used the same fitting procedures, but 
only allowed photoionization from the outer va- 
lence electron shell). 

The data set for collision strengths, transition 
probabilities and energy levels has been substan- 
tially updated using version 5.2 of the CHIANTI 
atomic database (Landi et al. 2006, Dere et al., 
1997). The updates are described in more detail 
in Section 12.31 Further updates are planned for 
the near future when version 6.0 of the database 
becomes available. 

Recombination lines from hydrogenic species are 
calculated using fits to the line emissivities of 
Storey & Hummer (1995). The data includes 
species up to Z=8, and temperatures extending 
to 10^ K for Z > 2. Species with Z > 8 are calcu- 
lated by scaling the He il data. The same hydro- 
genic scaling to He ii is also used for species with 
2 < Z < 8 when the local electron temperatures ex- 
ceed the highest tabulated temperature of 10^ K. 

Recombination lines from He-like species are cal- 
culated using fits to recombination data that will 
become publicly available in the next version 
of CHIANTI database. The available total (di- 
rect -|-cascade) temperature-dependent recombina- 
tion coefficients to the individual levels are used in 
the statistical equilibrium matrix that calculates 
the level populations, such that the final emissiv- 
ities automatically include contributions from re- 
combination. As well as the radiative-only rates 
derived from the Mewe et al. (1985) calculations 
for He-like iron, the dataset also includes more re- 
cent data of Porquet & Daubau (2000), which in- 



clude dielectronic contributions, when significant, 
for He-like ions with Z = 6, 7, 8, 10, 12 and 14. 

• Recombination contributions to Li-like Fe xxiv 
transitions are calculated in a similar fashion to 
He-like contributions (see item above), using the 
temperature-dependent recombination coefficients 
of Gu et al. (2003). The data, however, extend 
only to the 3d level, transitions originating from 
higher levels (e.g. Fe xxiv 7.99 A) are calculated 
via hydrogenic scaling of the coefficients given in 
Storey & Hummer (1995). 

• The ionization rate coefficients of atoms and ions 
by electron impact now uses the fits presented 
by Voronov (1997). The coUisional ionization 
and heating effects of the suprathermal secondary 
electrons following inner-shell photoionization are 
computed using approximations from ShuU & van 
Steenberg (1985). Ionization by suprathermal sec- 
ondaries is generally only important for very low 
ionization fractions. 

2.2.2. Inner-shell photoionization 

Previous versions of the MOCASSIN code were designed 
to treat the transfer of low-energy ionizing radiation, typ- 
ical of H II regions and PNe. In these environments, only 
ionization from the outer valence electron shell needs to 
be considered. We have now lifted this limitation by in- 
cluding inner shell photoionization using the Auger yields 
of Kaastra & Mewe (1993) and the photoionization cross- 
sections given by Verner et al. (1996). The solution of 
the ionization balance equations is rendered more compli- 
cated by the coupling of non- adjacent states through the 
emission of multiple electrons by high energy photons. 
We adopt a similar iteration scheme as that described 
by Ferland (2006) to solve the resulting ionization bal- 
ance matrix. 

2.2.3. Thomson/ Compton scattering 

The Compton cross-section is calculated via the Klein- 
Nishina formula, which simplifies to the Thompson case 
at non-relativistic energies. Compton heating and cool- 
ing contributions are calculated using standard formulae 
(e.g. Rybicki & Lightman, 1986). Bound Compton ion- 
ization and heating from high energy photons (^2.3 keV 
for hydrogen) are also taken into account, using the for- 
malism described by Ferland (2006, Eqn 294). 

The great advantage of a 3D Monte Carlo transfer 
technique is that scattering events can be treated taking 
into account the real geometry of the object. The re- 
processed directions of Compton-scattered packets, for 
example, are calculated stochastically using probability 
density functions based on the redistribution functions 
obtained from the Klein- Nishina formulae. This is of par- 
ticular interest for the calculation of fluorescence spectra, 
as described in Drake, Ercolano & Swartz (2007), Drake 
& Ercolano (2007a) and Drake & Ercolano (2007b), but 
it can in principle be applied to any line or continuum 
energy packet. 

2.3. Energy levels, collision strengths and radiative 
decay rates 
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Energy levels, collision strengths and radiative decay 
rates have been updated with data from version 5.2 of 
the CHIANTI atomic database (Landi et al., 2006; Dere 
et al., 1997). Since highly excited levels are generally not 
significantly populated in photoionized plasmas, not all 
data from CHIANTI were converted into the MOCASSIN 
format. Thus, for example, while the CHIANTI models 
for many of the boron-like ions consist of 125 levels, only 
data for the lowest 10 levels are currently used by MO- 
CASSIN. The number of levels used for each iso-electronic 
sequence are listed in Table [TJ The Maxwellian-averaged 
collisi on strengths, T, are stored as spline fits in CHI- 
ANTI (jPere et al.lll997( ) which allow T to be calculated 
for any temperature. For MOCASSIN, Ts are calculated 
from the CHIANTI spline fits for the temperature range 
2.0 < log T < 6.4 at 0.2 dex intervals. 

Up until version 4 of CHIANTI all collision data were 
fitted with a five point spline, following the method of 
I Burgess &: Tullvl (|1992l ). For some transitions, however, 
a five point spline was not able to accurately fit the en- 
tire set of collision strengths and thus it was necessary 
to omit some of the original data points for the fit. The 
choice of which data points to omit was influenced by 
the major application of CHIANTI - the analysis of emis- 
sion line spectra from collisionally-ionized plasmas such 
as the solar corona - thus care was taken to ensure the 
fits were accurate in the temperature range where the 
ion was most abun dant. (Tables of thes e temperature 
ranges are given in iMazzotta et anil998l for example.) 
Oftentimes this led to the collision strength data at low 
temperatures being omitted from the fit. With version 
4 of CHIANTI (jYoung et al.l I2003D nine point spline fits 
to the collision strength became available as an option, 
thus allowing improved fits at low temperatures. As it 
is the low temperatures that are important for photoion- 
ized plasmas, the early CHIANTI data-sets have been re- 
assessed to ensure that the low temperature data points 
are being accurately reproduced. Poor fits were improved 
by re-fitting the original data with nine point splines. An 
example of an improved fit to the O IV ground transition 
is shown in Fig. [TJ 

In the process of adapting the CHIANTI data files for 
MOCASSIN, the existing data within MOCASSIN have been 
compared with those in CHIANTI. In general, MOCASSIN 
used the same or older atomic data; in the latter case 
CHIANTI data simply replaced the MOCASSIN data (after 
the checks described in the previous paragraph had been 
performed). In some cases, MOCASSIN made use of bet- 
ter atomic data (either more recent calculations, or data 
more suited to photoionized plasmas), and so these data 
were assessed and added to CHIANTI. The comparison 
between CHIANTI and MOCASSIN also proved valuable for 
identifying errors in the data files and a number of cor- 
rections have been made to chianti. These updates will 
appear in the next version of that database. 

3. BENCHMARK TESTS 

The complexity of the calculations involved in large 
scale photoionization codes, where a number of coupled 
microphysical process are at play, imposes the need for 
careful testing before predictions from such codes may 
be trusted. In terms of energetics, care should be taken 
that all ionization and recombination, heating and cool- 
ing channels are correctly represented. 



The new code presented here has been thoroughly 
tested and compared, for the ID case, to the solutions 
to a number of benchmark models obtained by inde- 
pendent codes routinely used in the literature. The set 
of benchmarks used was presented by Pcquignot et al. 
(2001) and comprises a plane-parallel simulation for con- 
ditions typical of NLRs of AGNs and three low-density 
thin slab models irradiated by a strong broken power- 
law continuum radiation field, as described in Table 10 of 
Pequignot et al. (2001), and reproduced here for conve- 
nience (see Table[3]). For the NLR model the illuminating 
spectrum is a power law in energy units with slope 1.3, 
and the ionization parameter, J/is.e = 0.01 is defined as 
-^13. 6 / (c X rifj), where -Fia.e is the incident photon flux 
above 13.6 eV. The ionization parameters of the three 
X-ray slab models, XO.l, Xl.O and XIO.O are f/ioo = 0.1, 
1.0 and 10.0, respectively, where the definition of C/ioo is 
analogous to that of [/13.6 above. 

The plane-parallel geometry assumed by the ID test, 
can be mimicked with our 3D code by modelling a thin, 
elongated cuboidal grid with plane parallel diffuse illu- 
mination coming at normal incidence from one of the 
smaller sides; the energy packets are only allowed to es- 
cape from the side opposed to the illuminating side and 
all other sides act as mirrors (for a discussion of the mir- 
ror technique see Ercolano et al., 2003b). The main input 
parameters for the benchmark tests are given in Table 1 
of Pequignot et al. (2001) and repeated in our Table [2] 
for convenience. 

We run a number of tests to estimate the random er- 
ror due to Monte Carlo sampling in order to establish an 
efficient, yet reliable, combination of cell numbers and en- 
ergy packets needed to sample them. We adopted grids 
of 300 depth points and three by three in the short axes, 
with the number of packets needed to achieve conver- 
gence varying from one to eight million. Errors of the 
order of 0.5 % in computed line strengths are commmon 
with this setup, with larger errors (<6%) sometimes ob- 
tained for the less abundant ions. 

We note that radiative transfer effects are largely unim- 
portant for the lines listed in the X-ray slab models. Fur- 
thermore, the models are optically thin to ionizing radi- 
ation and electron temperature and density, and ioniza- 
tion structure variations through the slabs are minimal. 

Before discussing the benchmark results in detail, we 
note that many atomic data updates and code devel- 
opments have occurred in the past seven years (see Er- 
colano 2005). In particular, the radiative and dielec- 
tronic recombination coefficients calculated by Badnell 
et al. (2003), Badnell (2006a,b), Zatsarinny et al. (2003, 
2004a,b, 2006), Colgan et al. (2003, 2004), Altun et 
al. (2004, 2006), Mitnik & Badnell (2004), covering aU 
elements up to Zn (plus Kr, Mo and Xe), including se- 
quences up to Na-like electron target^, represent a sig- 
nificant improvement over earlier sets. 

3.1. Narrow Line Region 

We have organized the benchmark comparison in a 
similar fashion to Pequignot et al. (2001); for the NLR 
model we check that mocassin's prediction of each quan- 
tity falls within the range of the values obtained by the 

^ The data set is publicly available form 
|http: //amdpp. phys.strath.ac.uk/tamoc/DATA/ 1 
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other codes, hsted in Table 10 of Pequignot et aL (2001). 
In those cases where mocassin's prediction hes outside 
the range, we compute the isolation factors, defined as 
the ratio of the value given by MOCASSIN to the upper 
limit of the range, if our value is above the range, or the 
ratio of the lower limit of the range to the value given by 
MOCASSIN, if our value is below the range. Table S] lists 
line fluxes, temperature and ionized helium fraction as 
predicted by MOCASSIN (column 4) as well as the mini- 
mum and maximum values (columns 2 and 3) predicted 
by the other codes and given in Table 7 of Pequignot et 
al (2001). The isolation factors computed for MOCASSIN 
are given in the last column of our Table [H where we 
have added a minus or plus sign to indicate whether our 
prediction lies below or above the range predicted by the 
other codes. 

In Pequignot et al. (2001) it is stated that isolation 
factors larger that 1.3 should be considered "indicative 
of a very significant departure and possible problem" . 
From Table S] we find that only six out of 48 quantities 
predicted carried an isolation factor with absolute value 
larger that 1.3, with only one of these being larger that 
2.0. We judge this to be indicative of a high degree of 
consistency between the results of mocassin and other 
codes, especially in the light of the atomic physics up- 
dates that have occurred between the year 2001 and to- 
day. Furthermore, the temperature and ionization struc- 
ture of our model is also in good agreement with those 
calculated by other codes indicating that the differences 
in some line predictions can probably be ascribed to 
atomic physics changes. 

The largest isolation factors were obtained for sulphur 
ions. The values in Table 3] show that our predictions 
for the [S iii] lines are generally lower than those pre- 
dicted by other codes, and we suspect that dielectronic 
recombination effects may be at play in this case. Rates 
for third row elements are poorly known. In the case of 
recombination of S IV into S ill the rates calculated by 
Nahar (2000) are available. These calculations used an 
ab-initio method to derive total (electron-|-ion) recombi- 
nation (Nahar & Pradhan 1994, 1995) which enables the 
consideration of radiative and dielectronic recombination 
processes in a unified manner. However not all model- 
ers choose to include the total recombination coefficients 
due to the uncertainties inherent to the method that re- 
lies on theoretical predictions of the resonance positions. 
Here we have used the data of Colgan et al. (2003) , which 
were not available to the 2001 benchmark modelers. Fur- 
thermore, no data was available in 2001 for S ill ^ S ii 
dielectronic recombination and estimates were used in- 
stead (e.g. Ali et al., 1991). Here we used the data of 
Colgan et al. (2004). It seems likely that some of the 
scatter recorded for these lines is indeed due the differ- 
ent atomic data available and assumptions made. 

As a final note we should add that the temperature at 
the inner edge calculated by MOCASSIN is slightly lower 
than the values obtained by the other codes. This is to be 
expected since MOCASSIN uses a exact treatment of the 
radiative transfer, including the diffuse component, while 
comparison codes approximate the transfer of the diffuse 
component by either iterating along only one direction or 
by adopting an "outward-only" approximation, but along 
several directions. As noticed by Pequignot et al. (2001), 
the kinetic temperature tends to be lower in the inner- 



most layers of models with exact transfer, as the ionizing 
radiation field there is softer. Figure [2] shows electron 
temperatures, ^^^[iir], and densities, no[cm~'^], calculated 
with MOCASSIN as a function of column density in the 
slab. 

3.2. X-ray irradiated plane-parallel slabs 

The Pequignot et al. (2001) benchmarks constitute 
the first attempt to assess the accuracy of photoioniza- 
tion codes in the X-ray regime. As well as atomic data 
issues, already pointed out in the previous section, a fur- 
ther complication for this set of models is posed by the 
fact that, in general, for the far-UV lines listed, the orig- 
inal benchmarks considered sums of multiplet lines in a 
rather liberal sense, so that energetics were privileged at 
the expense of accuracy of wavelengths (with the excep- 
tions of some fine-structure optical lines). It is therefore 
possible that part of the scatter amongst the original 
benchmark values arose because of the different way mul- 
tiplets have been handled in different model calculations 
(Pequignot, priv. comm.). This problem persists and 
should be kept in mind when looking at the comparisons 
performed in this article. The original benchmark tests 
were quite rough and their main aim was to highlight 
gross discrepancies, particularly in derived temperatures. 
A new benchmarking exercise has yet to appear for the 
high-energy regime; however it is still useful for us to 
check our code against the 2001 benchmarks as major 
problems in the thermal and ionization balances would 
be uncovered by such a comparison. Furthermore, we 
propose that the solutions presented here be taken as an 
improved set, obtained with up-to-date atomic data. 

Models of X-ray slabs are assessed by comparing, for 
each quantity, mocassin's prediction with the minimum 
and maximum limits of the range of values obtained by 
the other codes and listed in Tables 11-13 of Pequignot 
et al. (2001). 

In Figures [Hill] and El we use red (striped) histograms to 
illustrate mocassin's results, while the filled and empty 
black histograms respectively mark the lower and the 
higher limit of the range of solutions obtained by the 
other codes in the original exercise. All line fiuxes are in 
units of [erg s~^] for a column density of 10^^ cm~^, with 
the emissivities being summed over 47r. In aid of future 
benchmarking efforts we provide in Table[5]the individual 
line fluxes in [erg s^^] as predicted by this version (3.00) 

of MOCASSIN. 

In spite of all the complications listed above, the agree- 
ment of mocassin's predictions with the benchmarks is 
reassuring and adds confidence to current and future as- 
trophysical applications of the code. In particular, very 
good agreement is shown for the relative strengths of 
lines from far infrared to X-ray wavelengths. Some ex- 
ceptions are discussed in more detail below. 

The agreement between MOCASSIN and the codes in 
Pequignot et al. (2001) for the benchmark XOl is good, 
although the electron temperature of the slab, listed in 
Table El carries an isolation factor of 1.14. While this 
is small and not of concern per se it may affect the 
strengths of some of the emission lines. 

Figures|4]and[5]demonstrate a very good agreement be- 
tween MOCASSIN and the codes in Pequignot et al. (2001) 
for most of the lines listed in Table [5] for benchmarks XI 
and XIO. The small discrepancies with some of the H-like 
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TABLE 1 

The number of energy levels used for each iso-electronic 
sequence when compiling the mocassin atomic data files. 



Iso-electronic 


Number 


Iso-electronic 


Number 


Sequence 


of levels 


Sequence 


of levels 


H 


9 


Na 


3 


He 


17 


Mg 


5 


Li 


15 


Al 


12 


Be 


10 


Si 


27 


B 


10 


P 


5 


C 


15 


S 


5 


N 


13 


CI 


2 


O 


9 


Ar 




F 


2 


K 


2 


Ne 




Ca 


9 






Sc 


19 



recombination lines (e.g. He ii 303.8 A and H i 1216 A 
in XI and Ar XVIII 20.20 A, N vil 24.78 A, C vi 28.47 A 
and 33.74 A and He ii 303.8 A in XIO) are most probably 
due to different extrapolation techniques of the Storey & 
Hummer (1995) data to the high temperatures of these 
models. 

4. SUMMARY 

We have presented a new version of the fully 3D Monte 
Carlo photoionization code, MOCASSIN. The code was 
extended to allow the modelling of plasma irradiated by 
a hard continuum spanning from radio to gamma rays. 
The atomic data set of the code was also significantly up- 
dated and it is now synchronized with the latest release 



of the CHIANTI database. The applicability and limita- 
tions of the new code were discussed, and the results of 
a thorough benchmarking exercise presented. 

No major problems were found by the benchmark tests, 
although some minor differences have been found. We 
have highlighted a number of significant improvements 
in the atomic datasets available today compared to those 
available at the time the original benchmark tables were 
compiled. We provide here updated values and we em- 
phasize the need of a new benchmarking exercise to be 
undertaken by the plasma modelling community. 

The good performance of the MOCASSIN code in all 
benchmark tests demonstrates that it is ready for appli- 
cation to real astrophysical problems. The public version 
of the X-ray enhanced MOCASSIN code is available on re- 
quest from the author. 
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Ne/H X 10^ 


10.0 
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11.0 
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Mg/H X 10^ 


3.0 
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3.7 


3.7 


Si/H X 10^ 


3.0 


3.5 


3.5 


3.5 


S/H X 105 


1.5 


1.6 


1.6 


1.6 


Ar/H X 10^ 




0.37 


0.37 


0.37 


Fe/H X 10^ 




4.0 


4.0 


4.0 
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TABLE 4 

Standard Narrow Line Region (NLR). ^ 



Quantity 


Min 


Max 


MOCASSIN 


isolation 


H/3 crg/s/cm'^ 


1.09 


1.49 


1.24 





H/3 4861 


1.00 


1.00 


1.00 





He I 5876 


0.100 


0.139 


0.111 





He II 4686 


0.226 


0.260 


0.226 


0. 


C II ]2325+ 


0.362 


0.600 


0.438 





CII1335 


0.09 


0.148 


- 





CIH] 1907+1909 


2.33 


4.09 


3.20 


0. 


CIV1549+ 


3.36 


4.90 


3.84 


0. 




NI]5200+5198 


0.034 


0.230 


0.033 


-1.03 




Nil] 6584+6548 


1.19 


3.27 


0.933 


-1.28 




Nil 5755 


0.018 


0.310 


0.013 


-1.39 


Nni]1749+ 


0.177 


0.265 


0.157 


-1.12 


[Nin]57.3/im 


0.042 


0.050 


0.041 


-1.03 


NIV] 1487+ 


0.203 


0.250 


0.343 


+1.37 


NV1240+ 


0.106 


0.191 


0.096 


-1.10 




OI]63.1/tm 


0.220 


3.10 


0.156 


-1.43 




OI] 6300+6363 


1.37 


2.17 


1.44 







Oil] 3726+3729 


1.25 


1.85 


1.23 


-1.02 




OIII]51.8/im 


0.667 


0.858 


0.919 


+1.07 




OIII]88.3/im 


0.092 


0.112 


0.117 


+1.04 




OIII] 5007+4959 


31.2 


34.9 


39.3 


+1.12 




OIII]4363 


0.286 


0.348 


0.273 


-1.04 




OIV]25.9Atm 


1.72 


2.53 


1.948 





OIV]1403+ 


0.311 


0.510 


0.374 





OV]1218+ 


0.120 


0.201 


0.148 





OVI1034+ 


0.030 


0.059 


0.041 







NeII]12.8/im 


0.145 


0.832 


0.165 







NeIII]15.5Atm 


1.38 


2.86 


1.44 







Nelll] 3869+68 


1.64 


2.46 


1.92 







NeIV]2423+ 


0.394 


0.586 


0.387 


-1.02 




NeV] 3426+3346 


0.520 


0.670 


0.456 


-1.14 




NeV] 24.2^4111 


0.162 


0.450 


0.340 







NeVI]7.63/im 


0.226 


0.288 


0.243 





MgH 2798+ 


1.23 


2.33 


1.77 


0. 




MgIV]4.49/im 


0.092 


0.116 


0.060 


-1.53 




MgV]5.6lAtm 


0.100 


0.230 


0.170 







SiII]34.8Atm 


0.728 


1.37 


0.664 


-1.09 


Sill] 2335+ 


0.114 


0.257 


0.177 





Siin]1892+ 


0.089 


0.312 


0.130 
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0.078 


0.121 


0.101 
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0.750 


1.86 


1.80 
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0.082 


0.180 


0.222 
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0.560 


0.968 


0.284 
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0.488 


0.140 
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2.32 


1.30 
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SIV]10.5Atm 


0.850 


1.74 


1.57 







■Tinner [-^1 


16840 


17100 


16250 


-1.04 


< T[N^+N,] > [K] 


11970 


12920 


12480 





<He+/He>/<H+/H> 


0.729 


0.766 


0.736 
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Fig. 1. — Plot demonstrating the difference between five and nine 
point spline fits for the Oiv ^P\/2 — ^^3/2 ground transition. The 
original collision strength data are from lZhang et al.l (l994). The 
five point spline fit was optimized to the temperatures relevant to 
an electron ionized plasma (~ 10^ K) and it fails to reproduce the 
low-temperature region. The nine point spline gives a good fit to 
the low temperature data points needed for photoionizcd plasmas. 
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Fig. 2. — NLR benchmark; electron temperature, Ta[K], and 
density, nc[cm~^], as a function of hydrogen column density, 
A^H [cm~^], through the slab. 

X01 




Fig. 3. — X-ray slab; Ux=0.1 (XOl). The red, striped histograms 
illustrate MOCASSIN'S results, while the filled and empty histograms 
represent the lower and higher limits of the range of predictions 
obtained by the codes in Pequignot et al. (2001) . 
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TABLE 5 

X-RAY IRRADIATED SLABS. EMISSION LINE FLUXES PREDICTED BY MOCASSIN, GIVEN IN UNITS OF [erg/s]. 





XOl 






XI 






XlO 




Ion 


A [A] 


Flux 


Ion 


A [A] 


Flux 


Ion 


A [A] 


Flux 


Mg9 


9.314 


1.17E-04 


S16 


4.792 


7.05E-05 


Fc26 


1.392 


1.9SE-05 


NclO 


12.13 


1.30E-04 


S15 


5.101 


1.40E-04 


Fc26 


1.425 


4.06E-05 


Nc9 


13.45 


1.29E-04 


Sill4 


6.182 


2.30E-04 


Fe26 


1.503 


l.llE-04 


Ne9 


13.55 


1.53E-04 


Sill3 


6.648 


7.32E-05 


Fe25 


1.573 


7.75E-05 


Ne9 


13.70 


4.27E-04 


Sill3 


6.688 


1.08E-04 


Fe26 


1.780 


6.23E-04 


08 


15.18 


1.19E-04 


Sill3 


6.740 


2.31E-04 


Fe25 


1.851 


5.72E-04 


08 


16.01 


3.27E-04 


Mgl2 


7.106 


4.91E-05 


Fe25 


1.859 


7.81E-04 


08 


18.97 


1.84E-03 


Mgl2 


8.421 


2.38E-04 


Fe25 


1.S6S 


1.29E-03 


N7 


20.91 


3.86E-05 


Mgll 


9.314 


8.98E-05 


Arl8 


3.151 


1.79E-05 


07 


21.60 


4.57E-04 


NelO 


10.24 


9.04E-05 


Arl8 


3.731 


8.75E-05 


07 


21.80 


4.78E-04 


NelO 


12.13 


4.14E-04 


816 


3.991 


5.76E-05 


07 


22.10 


1.58E-03 


Ne9 


13.70 


5.61E-05 


Arl7 


3.994 


2.23E-05 


N7 


24.78 


2.07E-04 


08 


15.18 


1.15E-04 


816 


4.729 


2.79E-04 


C6 


28.47 


9.70E-05 


08 


16.01 


2.94E-04 


8il4 


4.947 


3.17E-05 


N6 


29.54 


1.24E-04 


08 


18.97 


1.26E-03 


815 


5.101 


3.54E-05 


C6 


33.74 


5.04E-04 


N7 


24.78 


9.64E-05 


8il4 


5.217 


S.27E-05 


C5 


41.47 


8.10E-05 


C6 


33.74 


1.68E-04 


8il4 


6.182 


3.7SE-04 


OS 


102.5 


1.94E-04 


Fel9 


101.5 


3.23E-05 


8il3 


6.740 


2.34E-05 


C6 


182.2 


5.01E-05 


08 


102.5 


1.05E-04 


Mgl2 


7.106 


4.7SE-05 


He2 


256.0 


2.32E-04 


Fel9 


108.4 


1.28E-04 


Mgl2 


8.421 


2.10E-04 


Fel5 


284.2 


6.87E-04 


Fe22 


101+17+36 


4.33E-04 


Fe26 


9.652 


6.60E-05 


He2 


303.8 


9.38E-04 


Fe20 


121.8 


3.03E-04 


NelO 


10.24 


6.57E-05 


Sill 


303.9 


2.77E-04 


Fe21 


128.0 


6.14E-04 


Fe25 


10.32 


7.46E-07 


Fel3 


316.0 


8.33E-05 


Fe20 


132.8 


4.44E-04 


Fe24 


10.63 


1.05E-04 


Fel4 


333.7 


1.39E-04 


Fe23 


132.8 


3.64E-04 


Fe24 


11.17 


3.37E-04 


Fel6 


335.4 


1.52E-04 


Fe24 


192.0 


7.72E-05 


NelO 


12.13 


2.81E-04 


Fel6 


360.8 


9.27E-05 


Fe24 


255.1 


3.90E-05 


08 


16.01 


6.57E-05 


Mg9 


368.1 


5.32E-04 


He2 


303.8 


1.64E-04 


OS 


18.97 


2.81E-04 


Ne7 


465.0 


3.22E-04 


Arl6 


389+ 


9.34E-05 


ArlS 


20.20 


2.34E-05 


Sill2 


506+ 


3.05E-04 


S14 


418+ 


1.58E-04 


N7 


24.78 


5.77E-05 


Sill 


581.0 


8.21E-05 


Sil2 


506+ 


7.30E-05 


816 


25.58 


2.64E-05 


MglO 


615+23 


1.04E-03 


Fe20 


567.8 


4.82E-05 


06 


28.47 


2.41E-05 


05 


629.7 


6.30E-05 


Fe20 


721.4 


9.57E-05 


8il4 


33.42 


3.3SE-05 


Ne8 


770+80 


1.47E-03 


Fe22 


845.4 


8.52E-05 


06 


33.74 


9.59E-05 


Ne7 


895.2 


5.39E-05 


HI 


1216 


9.13E-05 


Mgl2 


45.51 


1.81E-05 


06 


1034 


8.74E-04 


Fe21 


1354 


2.14E-04 


NelO 


65.56 


2.34E-05 


HI 


1216 


4.97E-04 


T/IO^K 




6.99 


OS 


102.5 


5.89E-05 


Fel3 


2579 


6.12E-06 






Fe23 


132.8 


1.42E-04 


Fel3 


3388 


5.60E-06 








Fe24 


192.0 


3.09E-04 


Fel4 


5303 


9.21E-04 








Fe24 


255.1 


1.3SE-04 


S12 


7536 


1.02E-04 








Hc2 


303.8 


7.97E-05 


FelS 


10747 


7.91E-05 








T/IO^K 




1.44 


Fel3 


10798 


5.22E-05 












SilO 


14302 


2.95E-04 














T/IO^K 




12.2 

















Fig. 4. — X-ray slab; Ux=l (Xl). The red, striped histograms 
illustrate MOCASSIN'S results, while the filled and empty histograms 
represent the lower and higher limits of the range of predictions 
obtained by the codes in Pequignot et al. (2001). 



Fig. 5. — X-ray slab; Ux=10 (XlO). The red, striped histograms 
illustrate MOCASSIN'S results, while the filled and empty histograms 
represent the lower and higher limits of the range of predictions 
obtained by the codes in Pequignot et al. (2001) . 
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